###########################################################################
# Lipps & Sczepanski & Malet
# Understanding Preferences over Borders
# International Studies Quaterly

# Appendix: Test effect of border distance on Likert scale ----------------
###########################################################################

library(tidyverse)
library(ggplot2)
library(viridis)
library(texreg)
library(raster)
library(sp)
library(sf)
library(geosphere)

#Load data
dat_long <- read.csv("survey_dat_clean_conjoint.csv", na.strings=c("","NA"), row.names = 1)
dat <- read.csv("survey_dat_clean.csv", na.strings=c("","NA"), row.names = 1, fileEncoding = "UTF-8")

#Recoding to dummy and bins
dat <- dat %>% mutate(dist_dummy = as.factor(case_when(distance_closest_border<=50 ~ 1,
                                                       distance_closest_border>50 ~ 0)),
                      borderdist = as.factor(case_when(distance_closest_border<=50 ~ "<50km",
                                                       distance_closest_border>50 & distance_closest_border<=100 ~ "<100km",
                                                       distance_closest_border>100 & distance_closest_border<=150 ~ "<150km",
                                                       distance_closest_border>150 ~ ">150km")))

dat$borderdist <- ordered(dat$borderdist, levels = c("<50km", "<100km", "<150km", ">150km"))
dat$ironcurtaindist <- ordered(dat$distance_ironcurtain, levels = c("<50km", "<100km", "<150km", ">150km"))

#create dummies for vote intention
pty <- fastDummies::dummy_cols(dat$vote_intention)
dat <- cbind(dat, pty[,c(2,3,5,6,9,10)])

dat <- rename(dat,
              AfD=".data_AfD",
              CDU=".data_CDU/CSU",
              FDP=".data_FDP",
              Greens=".data_Greens",                
              SPD=".data_SPD",                 
              Left=".data_The Left") #the option "no answer" is dropped due to singularities

dat$age <- as.numeric(dat$age)

#recode long data frame
dat_long$border <- as.factor(dat_long$border)
dat_long$country <- as.factor(dat_long$country)
dat_long$aim <- as.factor(dat_long$aim)

dat_long <- dat_long %>% mutate(borderdist = as.factor(case_when(distance_closest_border<=50 ~ "<50km",
                                                                 distance_closest_border>50 & distance_closest_border<=100 ~ "<100km",
                                                                 distance_closest_border>100 & distance_closest_border<=150 ~ "<150km",
                                                                 distance_closest_border>150 ~ ">150km")))

#### TABLES AND FIGURES ####--------------------------------------------------------

#Figure E1: Distance to nearest border

pdf(file = "fig_E1.pdf",width = 8,height = 5,onefile=F)

ggplot(dat, aes(distance_closest_border)) +geom_histogram(binwidth = 10) +
  labs(x="Distance to nearest border in km", y="Frequency") +
  theme_bw() +
  theme(
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size=12),
    axis.title.y= element_text(size = 14),
    axis.title.x = element_text(size=14))

dev.off()

#Figure E2: Map

shape <- st_read("gadm41_DEU_2.shp")

shape$CC_2 <- as.integer(shape$CC_2)

sub <- dat %>% 
  group_by(krs_name) %>% 
  mutate(means_krs_open = mean(border_open_close_1, na.rm=TRUE),
         n_krs_respon = n()) %>%
  dplyr::select(c(means_krs_open, krs_code, n_krs_respon)) %>% 
  group_by(krs_code) %>% 
  unique() %>%
  ungroup()

sub <- sub[complete.cases(sub),]
shape <- left_join(shape, sub, by = c("CC_2"="krs_code"))

pdf(file = "fig_E2.pdf",width = 10,height = 16,onefile=F)

ggplot() +
  geom_sf(data = shape, aes(fill=means_krs_open) ) +
  scale_fill_continuous(breaks=c(0,9), labels=c("open", "close"), type="viridis", direction=-1,name="Border should be    ") +
  theme_minimal() +
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=14),
        legend.title = element_text(size=14))

dev.off()

summary(sub$n_krs_respon)

#Table E1: Likert-scale
# use a dummy for border distance
mod3 <- lm(border_open_close_1 ~ AfD+FDP+CDU+SPD+Left+Greens+
             #left_right_1 + immigration_general_1 + european_integration_1 +
             gender + age + urban_rural + educat + migration_background + dist_dummy,
           data=dat)
mod4 <- lm(border_open_close_1 ~ #voteint2 +
             left_right_1 + immigration_general_1 + european_integration_1 +
             gender + age + urban_rural + educat + migration_background + dist_dummy,
           data=dat)

mod5 <- lm(border_open_close_1 ~ AfD+FDP+CDU+SPD+Left+Greens+
             #left_right_1 + immigration_general_1 + european_integration_1 +
             gender + age + urban_rural + educat + migration_background + borderdist,
           data=dat)

mod6 <- lm(border_open_close_1 ~ #voteint2 +
             left_right_1 + immigration_general_1 + european_integration_1 +
             gender + age + urban_rural + educat + migration_background + borderdist,
           data=dat)

texreg(list(mod5,mod6),
       file = "tab_E1.tex",
       include.ci = F,
       digits = 3,
       single.row = TRUE,
       custom.header = list("Preference for closed border (0-10)" = 1:2))

#Figure E3: Conjoint by distance
results_border <- cregg::cj(chosen ~ border+aim+country, data=dat_long,
                            id=~id, estimate="mm", h0=0.5, by=~borderdist)

labels<- as_labeller(c("border" = "Measure", "aim"="Aim", "country"="Country", "sim_cat" = "Country"))

pdf(file = "fig_E3.pdf",width = 8,height = 5,onefile=F)

ggplot(results_border)+
  geom_point(aes(x=estimate, y=level, color=borderdist), size=2, position = position_dodge(0.4))+
  geom_errorbar(aes(xmin=lower, xmax=upper, y=level, color=borderdist), linewidth=1, position = position_dodge(0.4))+
  facet_grid(feature~., scales = "free", labeller = labels)+
  geom_vline(xintercept = 0.5, linetype="dotted")+
  labs(x="Marginal Means", y="", color="Closest land border")+
  scale_color_viridis(discrete=T)+
  theme_minimal() +
  theme(
    axis.text.y.left = element_text(size = 16),
    axis.text.x = element_text(size=16),
    axis.title.y= element_text(size = 16),
    axis.title.x = element_text(size=16),
    strip.text.y.right = element_text(size=16),
    legend.text = element_text(size=16),
    legend.position = "top")

dev.off()
